# Table B.5
rm(list=ls())

(WD <- getwd())
if (!is.null(WD)) setwd(WD)

indata0      = paste0(WD,"/_individual_data/_spfrawdata/","", collapse = NULL)

library(haven)
library(AER)
library(sandwich)
library(lmtest)
library(pracma)
library(stargazer)
library(plm)
library(pracma)
library(DataCombine)
library(jtools)
library(plyr)
library(ggplot2)
library(tidyverse)
library(dotwhisker)
library(DescTools)

nobs_win = 20

# PGDP
load(paste(indata0,"us_spf_pgdp_ind.Rda",sep=""))
data$ID         = data$id
data            = data[data$qdate>=1970.00,]
data            = data[data$qdate<=2020.00,]
data            = ddply(data,.(qdate),transform, cons = mean(f2,na.rm = TRUE))
data            = ddply(data,.(qdate),transform, fc_rev_avr = mean(fc_rev,na.rm = TRUE))

tmp             =  max(data$ID)-min(data$ID)+1
beta_ind_save   =  rep(NA, tmp)
iter            =  0
for (ii in min(data$ID):max(data$ID)){
  iter     = iter + 1
  data_tmp = data[data$ID==ii,]
  nobs     = nrow(data_tmp)
  if (nobs>= nobs_win) {
     tryCatch({plm_tmp    = plm(fc_err ~ fc_rev, data=data_tmp,index=c("ID", "qdate"), model="pooling")
     beta                 = coefficients(plm_tmp)
     beta_ind_save[iter]  = as.numeric(beta[1])}, error = function(e){})
      }
}
median_beta    = median(beta_ind_save, na.rm=TRUE)

tmp             =  max(data$ID)-min(data$ID)+1
delta_ind_save  =  rep(NA, tmp)
b_ind_save      =  rep(NA, tmp)
iter            =  0
for (ii in min(data$ID):max(data$ID)){
  iter     = iter + 1
  data_tmp = data[data$ID==ii,]
  nobs     = nrow(data_tmp)
  if (nobs>= nobs_win) {
    tryCatch({lm_tmp     = lm(fc_err ~ cons, data=data_tmp)
    delta                = coefficients(lm_tmp)
    delta_ind_save[iter] = as.numeric(delta[2])}, error = function(e){})
    
    tryCatch({lm_tmp     = lm(fc_err ~ fc_rev_avr, data=data_tmp)
    b                    = coefficients(lm_tmp)
    b_ind_save[iter]     = as.numeric(b[2])}, error = function(e){})
  }
}
median_delta    = median(delta_ind_save, na.rm=TRUE)
median_b        = median(b_ind_save, na.rm=TRUE)

median_est =  data.frame(median_b, median_beta, median_delta)
rownames(median_est) = c("Estimates")
median_est = t(median_est)
row.names(median_est) = c('Avr. Forecast Revision', 'Ind. Forecast Revision', 'Prev. Consensus')
write.table(format(median_est,digits=3), file = "_tables/table_b5_top.txt", sep = "\t")

# CPI
load(paste(indata0,"us_spf_cpi_ind.Rda",sep=""))
data$ID         = data$id
data            = data[data$qdate>=1970.00,]
data            = data[data$qdate<=2020.00,]
data            = ddply(data,.(qdate),transform, cons = mean(f2,na.rm = TRUE))
data            = ddply(data,.(qdate),transform, fc_rev_avr = mean(fc_rev,na.rm = TRUE))

tmp             =  max(data$ID)-min(data$ID)+1
beta_ind_save   =  rep(NA, tmp)
iter            =  0
for (ii in min(data$ID):max(data$ID)){
  iter     = iter + 1
  data_tmp = data[data$ID==ii,]
  nobs     = nrow(data_tmp)
  if (nobs>= nobs_win) {
    tryCatch({plm_tmp    = plm(fc_err ~ fc_rev, data=data_tmp,index=c("ID", "qdate"), model="pooling")
    beta                 = coefficients(plm_tmp)
    beta_ind_save[iter]  = as.numeric(beta[1])}, error = function(e){})
  }
}
median_beta    = median(beta_ind_save, na.rm=TRUE)

tmp             =  max(data$ID)-min(data$ID)+1
delta_ind_save  =  rep(NA, tmp)
b_ind_save      =  rep(NA, tmp)
iter            =  0
for (ii in min(data$ID):max(data$ID)){
  iter     = iter + 1
  data_tmp = data[data$ID==ii,]
  nobs     = nrow(data_tmp)
  if (nobs>= nobs_win) {
    tryCatch({lm_tmp     = lm(fc_err ~ cons, data=data_tmp)
    delta                = coefficients(lm_tmp)
    delta_ind_save[iter] = as.numeric(delta[2])}, error = function(e){})
    
    tryCatch({lm_tmp     = lm(fc_err ~ fc_rev_avr, data=data_tmp)
    b                    = coefficients(lm_tmp)
    b_ind_save[iter]     = as.numeric(b[2])}, error = function(e){})
  }
}
median_delta    = median(delta_ind_save, na.rm=TRUE)
median_b        = median(b_ind_save, na.rm=TRUE)

median_est =  data.frame(median_b, median_beta, median_delta)
rownames(median_est) = c("Estimates")
median_est = t(median_est)
row.names(median_est) = c('Avr. Forecast Revision', 'Ind. Forecast Revision', 'Prev. Consensus')
write.table(format(median_est,digits=3), file = "_tables/table_b5_bottom.txt", sep = "\t")
